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Abstract We propose a new approximate Bayesian computation (ABC) algorithm that aims at minimizing 
the number of model runs for reaching a given quality of the posterior approximation. This algorithm 
automatically determines its sequence of tolerance levels and makes use of an easily interpretable stopping 
criterion. Moreover, it avoids the problem of particle duplication found when using a MCMC kernel. When 
applied to a toy example and to a complex social model, our algorithm is 2 to 8 times faster than the three 
main sequential ABC algorithms currently available. 

Keywords Approximate Bayesian computation • Population Monte Carlo • Sequential Monte Carlo • 
Complex model 



1 Introduction 

Approximate Bayesian computation (ABC) techniques appear particularly relevant for calibrating stochas- 
tic models because their very principle includes stochasticity but they are applicable to any model. They 
generate a sample of model parameter values (0,),=i...,a? (often also called particles) from the prior distri- 
bution 7t{9) and select the 0, values leading to model outputs x ^ satisfying a proximity criterion 
with the target datay {p{x,y) < e, p expressing a distance, e being a tolerance level). The selected sample 
of parameter values approximates the posterior distribution of parameters, leading to model outputs with 
the expected quality of approximation. However, in practise, running these techniques is very demanding 
computationally because sampling the whole space of parameters requires a number of simulations which 
grows exponentially with the number of parameters to identify. This tends to limit the application of these 
techniques to easily computable models (Beaumont 2010). In this paper, our goal is minimizing the num- 
ber of model runs for reaching a given quality of posterior approximation, and thus to make the approach 
applicable to a larger set of models. 

ABC is the subject of intense scientific researches and several improved versions of the original scheme 
are available, such as using local regressions to improve parameter inference (Beaumont et al. 2002; Blum 
and Francois 2010), automatically selecting informative summary statistics (Joyce and Marjoram 2008; 
Fearnhead and Prangle 2011), coupling to Markov chain Monte Carlo (Marjoram et al. 2003; Wegmann 
et al. 2009) or improving sequentially the posterior distributions with sequential Monte Carlo methods 
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(Sisson et al. 2007; Toni et al. 2009; Beaumont et al. 2009). This last class of methods approximates pro- 
gressively the posterior, using sequential samples 5^'^ = (0/'^),=i,..,a? derived from sample S^'^^\ and using 
a decreasing set of tolerance levels {ei, ...,£7-}. This strategy focuses the sampling effort in parts of the 
parameter space of high likehhood, avoiding to spend much computing time in systematically sampling the 
whole parameter space. 

The first sequential method apphed to ABC was proposed by Sisson et al. (2007) with the ABC-PRC 
(Partial Rejection Control). This method is based on a theoritical work of Del Moral et al. (2006) to ABC. 
However, Beaumont et al. (2009) has shown that this method leads to a bias in the approximation of the 
posterior. Beaumont et al. (2009); Toni et al. (2009) proposed a new algorithm, called Population Monte 
Carlo ABC in Beaumont et al. (2009) and hereafter called PMC. This algorithm, corrects the bias by affect- 
ing to each particle a weight corresponding to the inverse of its importance in the sample. It is particularly 
interesting in our perspective because it provides with a rigorous framework to the sequential sample idea, 
which seems a good way for minimizing the number of runs. In this approach, the problem is then defining 
the sequence of tolerance levels {ci, ...,£7-}. Drovandi and Pettitt (201 1) and Del Moral et al. (2012) solve 
partly this problem by deriving the tolerance level at a given step from values of p{x,y) of the previously 
selected sample. However, a difficulty remains: when to stop? If the final tolerance level Et is too large, 
the final posterior will be of bad quality. Inversely, a too small Et leads to a posterior that could have been 
obtained with less model runs. 

Moreover, the MCMC kernel used in Drovandi and Pettitt (2011) and Del Moral et al. (2012) to sample 
new values Oj'\ despite its mathematical elegance, has a significant drawback in our view: it can lead 
to particle duplications. Indeed, each time the MCMC jumps from a particle to a new one which is not 
accepted, the initial particle is kept in the new sample of particles, hence when this occurs several times with 
the same initial particle, this particle appears several times in the new sample. As long as these duplications 
are very few compared with the size of the sample, their effect can be neglected, but it can easily happen 
that their number grows to a very significant part of the sample, then strongly deteriorating the quality of 
the posterior, as illustrated below. To solve this problem, Drovandi and Pettitt (2011) proposed to perform 
R MCMC jump trials instead of one, while Del Moral et al. (2012) proposed to resample the parameter 
values when too many are duplicated. Del Moral et al. (2012) also proposed to run the model M times for 
each particle, in order to decrease the variance of the acceptance ratio of the MCMC jump. However, these 
solutions increase the number of model runs, going against the initial benefit of using sequential samples. 

In this paper, we propose a modification of the PMC algorithm that we call adaptive population Monte 
Carlo ABC (hereafter called APMC). This new algorithm determines by itself the sequence of tolerance 
levels as in Drovandi and Pettitt (201 1) and Del Moral et al. (2012), and it also provides a stopping criterion. 
Furthermore, our approach avoids the problem of duplications. We prove that the computation of the weights 
associated to the particles in this algorithm lead to the intended posterior distribution and we also prove that 
the algorithm stops whatever the chosen value of the stopping parameter We show that our algorithm, 
applied to a toy example and to an individual-based social model, requires significantly less simulations to 
reach a given quality level of the posterior distribution than the PMC algorithm of Beaumont et al. (2009), 
the replenishment SMC ABC algorithm of Drovandi and Pettitt (2011) (hereafter called RSMC) and the 
adaptive SMC ABC algorithm of Del Moral et al. (2012) (hereafter called SMC). These algorithms are 
detailed in Appendix A. 



2 Adaptive population Monte-Carlo approximate Bayesian computation 

2.1 Overview of the APMC algorithm 

The APMC algorithm follows the main principles of the sequential ABC, and defines on-line the tolerance 
level at each step like in Drovandi and Pettitt (2011) and Del Moral et al. (2012). For each tolerance level 
it generates a sample 5^'^ of particles and computes their associated weights. This weighted sample 
approximates the posterior distribution, with an increasing approximation quality as £, decreases. We say 
that a parameter value 9^, satisfies the tolerance level £,, if when running the model we get x ~ f{x\9P), 
such that its distance pj'^ — p {x,y) to the target data y, is below £,. Suppose the APMC reached step f — 1, 

with a sample S^'^^'' of Na = [(XN\ particles and their associated weights {9j' '\wf '')/=i,..,A'a, the main 
features of the APMC are (see Algorithm 5 for details): 
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- the algorithm generates A^ — A^ce particles ' )j=Na+i..;N where df ' ' - ^ ( 0* , (T(^,_ i ) ) , the seed 9* 
is randomly drawn from the weig htedset(0/' ^\wf ''),=i,..,Ar„ and the variance ar, of the Gaus- 
sian kernel (J^^^ jj) is twice the empirical variance of the weighted set (s/' ^\wf '■'),=i,...Afo,, 
following Beaumont et al. (2009). 

- the weights of the new particles {9j' ^^)j=Na+i,--,N are computed so that these new particles 
can be combined with the sample of the previous step without causing a bias in the posterior 
distribution. These weights are given by Eq. 2 (see below). 

- the algorithm concatenates the A^^ previous particles {9^' ^^)i=i Na '^i^h the N — Na new particles 

j=Na+i,...,N, together with their associated weights and distances to the data. This constitutes a 
new set noted S^ll,p = {9^ ,wf\pP)i=i,.,.N. 

- the next tolerance level £, is determined as the first a— quantile of the (p/'^),=i,..,a?. 

- the new sample 5^ = (0/'\wf'),=i,. jv„ is then constituted from the A^^ particles of Sf2np satisfying the 
tolerance level e,. 

- if the proportion pacc of particles satisfying the tolerance level i among the A^ — Na newly generated 
particles is below a chosen value Pacc^^,,, the algorithm stops, and its result is {9j''^)i=i^,,^Na with their 
associated weights. 

Note that in our algorithm, to get a number Na of retained particles for the next step, the choice of e, 
is heavily constrained: it has to be at least equal to the first a— quantile of the {pP)i=i,...N and smaller 

than the immediately superior (pP) value. We chose to fix it to the first a— quantile for simplicity. This 
choice also ensures that the tolerance level decreases from one iteration to the next: in the worst case where 
Pace = (no newly simulated particles accepted), e, = e,_ i . Our algorithm does not use a MCMC kernel 
and avoids duplicating particles. It requires a reweighting step in 0{Na) instead of 0{Na) in Drovandi and 
Petti tt (201 1), but in our perspective, this computational cost is supposed negligible compared with the cost 
of running the model. 



2.2 Weights correcting the kernel sampling bias 

As pointed out by Beaumont et al. (2009), the newly generated particles 9^ in a sequential procedure are 

no more drawn from the prior distribution but from a specific probability density dP that depends on the 
particles selected at the previous step and on the chosen kernel. This introduces a bias in the procedure. This 
bias should be corrected by attributing a weight equal to n{9P)/dP to each newly generated particle 9f \ 
The density of probability to generate particle 9^ at step t is given by the sum of the probabilities 
to reach 9^ from one of the Na particles of the previous step times their respective weights: 

where (p{x) = -p^e i is the kernel function. 

This yields the expression of the weight to be attributed to the newly drawn particle 9^: 

J') ^ ^ n) 

' "i;i^.(wr^Vlff,vvr'')a-\.p(a-\(e<')-ej'-^)) 

This formula differs from the scheme of Beaumont et al. (2009) where the weights need only to be 
proportional to Eq. 2 at each step. Since we want to concatenate particles obtained at different steps of the 
algorithm (while Beaumont et al. (2009) generate the sample at step t from scratch), we need the scaling of 
weights to be consistent across the different steps of the algorithm. Using the weight of Eq. 2 guarantees the 
correction of the sampling bias throughout the APMC procedure and ensures that the Na weighted particles 
9^ produced at the f-th iteration follow the posterior distribution n{9\p{x,y) < Et). 
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2.3 The stopping criterion 

We stop the algorithm when the proportion of "accepted" particles (Eq. 3) among the — Na new particles 
is below a predetermined threshold Pucc,,,,,, ■ This choice of stopping rule ensures that additional simulations 
would only marginally change the posterior distribution. Note that this stopping criterion will be achieved 
even if Pacc„,i„ = 0, this ensures that the algorithm converges. We present a formal proof of this assertion in 
Appendix B. 

3 Experiments on a toy example 

We consider four algorithms: APMC, PMC, the SMC and the RSMC. Their implementations in R (R De- 
velopment Core Team 2011) are available We compare them on the toy example studied in Sisson et al. 
(2007) where 7r(0) = 'i'[_io_io] and/(jt;|0) - (0, j^) + ^0 (0, 1) where (^u,^^) is the normal density 
of mean and variance O^. In this example, we consider thaty = is observed, so that the posterior density 
of interest is proportional to (0 (O, j^) + (j) (0, 1)) 7i{6). 

We structure the comparisons on two indicators: the number of simulations performed during the ap- 
plication of the algorithms, and the L2 distance between the exact posterior density and the histogram of 
particle values obtained with the algorithms. This L2 distance is computed on the 300-tuple obtained by 
dividing the support [—10, 10] into 300 equally-sized bins. 

We choose = 5000 particles and a target tolerance level equal to 0.01. For the PMC algorithm we use 
a decreasing sequence of tolerance levels from £\ =2 down to eu = 0.01. For the SMC algorithm, we use 
3 different values for a: {0.9,0.95,0.99} andM = 1 as in Del Moral et al. (2012). For the RSMC algorithm 
we use a = 0.5 as in Drovandi and Pettitt (2011). To explore our algorithm, we test 9 different values 
for a: {0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9}, and 4 different values for Pacc„„„- {0.01,0.05,0.1,0.2}. In 
each case, we perform 50 times the algorithm, and compute the average and standard deviation of the two 
indicators: the total number of simulations and the L2 distance between the exact posterior density and the 
histogram of particle values. We used as kernel transition a normal distribution parameterized with twice 
the weighted variance of the previous sample, as in Beaumont et al. (2009). 

We report below the effects of varying a and Pacc„,i„ on the performance of our algorithm, and compare 
it with the PMC, SMC and RSMC algorithms. 

3.1 Particle duplication in SMC and RSMC 

The number of distinct particles decreases during the course of the SMC algorithm whatever the value of 
a, as shown on Fig. la-b. The oscillations of the number of distinct particles are caused by the resampling 
step in the SMC algorithm (see Del Moral et al. (2012)), but they are not sufficient to counterbalance the 
overall decrease. This decrease deteriorates the posterior approximation as shown on Fig. 2. For the RSMC 
algorithm, the number of distinct particles is maintained at a reasonably high level (Fig. Ic), but this has 
a cost in terms of the number of required model runs (see Fig. 2). Note that the APMC and the PMC 
algorithms keep distinct particles. 

3.2 Influence of parameters on APMC 

The values of a and Pacc^j,, have an impact on the studied indicators. We find that smaller a and Pacc„,i„ 
improve the quality of the approximation (smaller L2 distance), and increase the total number of model runs, 
with Pacc„,i„ having the largest effect (Fig. 2). With a large a, the tolerance levels decrease slowly and there 
are numerous steps before the algorithm stops. In this toy example, our simulations show that all explored 
sets of (a , Pacc,„i„) such that Pacc,„i„ < 0.1 give good results for the criterion Number of simulations x 
(Fig. 3b). Large a provide slightly better results for small Pace,,,/,, while small a provide slightly better 



http : //motive . cemagref . f r/people/maxime . leiiormand/script_r_toyex 
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(a) (b) (c) 




1 2 3 4 5 10 15 20 25 30 35 5 10 15 20 25 30 35 



Number of simulations (xlO^) Number of simulations (xlO^) Number of simuiations (xlO^) 

Fig. 1 Number of distinct particles in a sample of N = 5000 particles during the course of the SMC and RSMC algorithms applied to 
the toy example, (a) SMC with a = 0.9 and M = 1; (b) SMC with a = 0.99 and M = 1; (c) RSMC with a = 0.5. In all three panels, 
the tolerance target is equal to 0.001. 



results for large Pacc„i„ (Fig- 3b). On this toy example it appears that intermediate values of a and Pacc„,i„ 
(0.3 < a < 0.7 and 0.01 < Pacc„,i„ < 0.05), present a good compromise between number of model runs and 
the quality of the posterior approximation. 
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Fig. 2 Posterior quality (L2) versus computing cost (number of simulations) averaged over 50 replicates. Vertical 
and horizontal bars represent the standard deviations among replicates. Algoiithm parameters used for APMC: a in 
{0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9} and Pacc„„ in {0.01,0.05,0.1,0.2}. Blue circles are used for Pacc„i^ =0.01, orange trian- 
gles for Pacc„,i„ = 0.05, green squares for Pace,,,,,, = 0. 1, and purple diamonds for Pacc,„i„ = 0.2. PMC: red plain triangles for a sequence 
of tolerance levels from ei = 2 down to en = 0.01. SMC: grey plain square for a in {0.9,0.95,0.99} (from left to right), M = 1 
and a e target equal to 0.01. RSMC: brown plain diamond for a = 0.5 and a e target equal to 0.01. Results obtained with a standard 
rejection-based ABC algorithm are depicted with black plain circles. 



3.3 Comparing performances 

Whatever the value of a and Pncc,,,/,,, 'he APMC algorithm always yields better results than the other 
three algorithms. It requires between 2 and 8 times less simulations to reach a given posterior quality 
L2 (Fig. 2). Furthermore, good approximate posterior distributions are very quickly obtained (Fig. 2). The 
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compromise between simulation speed and convergence level can also be illustrated using the criterion 
Number of simulations x Lj (Glynn and Whitt 1992). This criterion is smaller for the APMC algorithm 
(Fig. 3a). 




PMC SMC RSMC ABC APMC 

a 

Fig. 3 (a) Boxplot of the criterion "squared L2 distance times tlie number of simulations" for the different ABC algorithms. APMC: 
for a in {0.1,0.2,0.3,0.4,0.5,0.6,0.7.0.8,0.9} and po„„„„ = 0.01; SMC: for a in {0.9,0.95,0.99}, M = 1 and a e target equal to 
0.01 ; RSMC: for a = 0.5 and a e target equal to 0.01 ; ABC: for a e target equal to 0.01 ; PMC: for a sequence of tolerance levels from 
ei = 2 to Ell = 0.01. (b) Criterion "squared L2 distance times the number of simulations" in the APMC algorithm for the different 
values of a and Pacc,„„, ■ Each cell depicts the average of the criterion over the 50 performed replicates of the APMC. 



4 Application to the model Sim Villages 

In this section, we check if our algorithm still performs better than the PMC, the RSMC and the SMC when 
applied to an individual-based social model developed during the European project PRIMA^. The aim of 
the model is to simulate the effect of a scenario of job creation (or destruction) on the evolution of the 
population and activities in a network of municipalities. 



4.1 Model and data 

The model simulates the dynamics of virtual individuals living in 7 interconnected villages in a rural area of 
Auvergne (a region of Central France). A single run of the model SimVillages with seven rural municipali- 
ties takes about 1.4 seconds on a desktop machine (PC Intel 2.83 GHz). The dynamics include demographic 
change (aging, marriage, divorce, births and deaths), activity change (change of jobs, unemployment, in- 
activity, retirement), and movings from one municipality to another or outside of the set. The model also 
includes a dynamics of creation / destruction of jobs of proximity services, derived from the size of the local 
population. More details on the model can be found in Huet and Deffuant (2011). The individuals (about 
3000) are initially generated using the 1990 census data of the National Institute of Statistics and Economic 
Studies (INSEE), some of them are given a job type and a location for this job (in a municipality of the 
set or outside), they are organised in households living in a municipality of the set. The model dynamics is 
mostly data driven, but four parameters cannot be directly derived from the available data. They are noted 
Op for I < p < 4, described in Table 1 . 

We use our algorithm to identify the distribution of the four parameters for which the simulations, 
initialized with the 1990 census data, satisfy matching criteria with the data of the 1999 and 2006 census. 
The set of summary statistics {5,„}i<,„<m and the associated discrepancy measure used p„j are described 



^ PRototypical policy Impacts on Multifunctional Activities in rural municipalities - EU 7th Framework Research Programme; 
2008-201 1; https : //prima. cemagref . f r/the-project 
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in Table 2. We note S,„ the simulated summary statistics and S'^,, the observed statistics. The eight summary 
statistics are normalized (variance equalization) and they are combined using the infinity norm (Eq. 4): 

) (4) 

l<m<M 

We first generate a sample of length from the prior where [a, b] is available for each parameter in 

Table 1 , with a Latin hypercube (Camell 2009) and we select the best Na particles . To move the particles, we 
use as kernel transition a multivariate normal distribution parameterized with twice the weighted variance- 
covariance matrix of the previous sample (Filippi et al. 201 1). 

As in the section 3, we perform a parameter study and compare APMC with its three competitors. For 
APMC, a varies in ({0.3,0.5,0.7}) and pacv,„,„ in ({0.01,0.05,0.1,0.2}), and we set Na ^ 5000 particles. 
For the PMC, SMC and RSMC we also set — 5000 particles and a tolerance level target equal to 1 .4. 
The tolerance value e = 1 .4 corresponds to the average final tolerance value we obtain with APMC for 
Pacc„,j„ — 0.01. Note that otherwise this final tolerance is difficult to set properly and a worse choice for 
this value would have lead to worse performances of these algorithms. For the PMC algorithm, we use 
the decreasing sequence of tolerance levels {3,2.5,2, 1.7, 1.4}. For the SMC algorithm, we use 3 different 
values for the couple (a,M): {(0.9, 1), (0.99, 1) , (0.9, 15)}. For the RSMC algorithm we use a = 0.5, as 
in Drovandi and Pettitt (201 1). For each algorithm and parameter setting, we perform 5 replicates. 

We approximated posterior density (unknown in this case) with the original rejection-based ABC algo- 
rithm, starting with = 10, 000, 000, selecting 7890 particles below the tolerance level £ = 1 .4. 

To compute the L2 distance between posterior densities, we divided each parameter support into 4 
equally sized bins, leading to a grid of 4'^ = 256 cells, and we computed on this grid the sum of the squared 
differences between histogram values. 



Table 1 SimVillages parameter descriptions 



Parameters 


Description 


Range 


01 


Average number of children per woman 


[0,4] 


02 


Probability to accept a new residence for a household 


[0,1] 


03 


Probability to make couple for two individuals 


[0,1] 


04 


Probability to split for a couple in a year 


[0,0.5] 



Table 2 Summary statistic descriptions 



Summary statistic 


Description 


Measure of discrepancy 


Si 


Number of inhabitants in 1999 


L| distance 


S2 


Age distribution in 1999 


X' distance 


S3 


Household type distribution in 1999 


distance 


S4 


Net migration in 1999 


Li distance 


Ss 


Number of inhabitants in 2006 


L) distance 


Si 


Age distribution in 2006 


X' distance 


Si 


Household type distribution in 2006 


X^ distance 


Ss 


Net migration in 2006 


Li distance 



4.2 Study of APMC result 

APMC yields a unimodal approximate posterior distribution for the model SimVillages (Fig. 4). Interest- 
ingly, parameters 9i and 64 are slightly correlated (Fig. 4c). This is logical since they have contradictory 
effects on the number of child in the population. What is less straightforward is that we are able to partly 
tease apart these two effects with the available census data, since we get a peak in the approximate posterior 
distribution instead of a ridge. 
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0.2 0.4 0.6 0.8 0.2 0.4 0.6 0.8 0.2 0.4 0.6 0.8 



92 62 63 

Fig. 4 Contour plot of the bivariate joint densities of 6, and 8j obtained with our algorithm, and with a = 0.5 and Pacc,„i„ = 0.01 ; (a) 
01 and 62; (b) 9, and 63; (c) 61 and 64; (d) 01 and 63; (e) 62 and 64; (f) 63 and 64. 

4.3 Influence of parameters on APMC 

As for the toy example, we find that the intermediate values of {a,pacc,„j„) that we used lead to similar 
results (Fig. 5c). In practice, we therefore recommend to use a = 0.5 and Pacc,„j„ between 0.01 and 0.05 
depending on the wished level of convergence. 



4.4 Comparing performances 

APMC requires between 2 and 7 times less simulations to reach a given posterior quaUty than the other algo- 
rithms L2 (Fig. 5a). Again, the gain in simulation number is progressive during the course of the algorithm. 
The Number of simulations x criterion is again smaller for the APMC algorithm (Fig. 5b). 



0.20 - 



0.05 - 




Number of simulations (xlO^) PMC SMC RSMC APMC 



Fig. 5 (a) Posterior quality (L2) versus computing cost (number of simulations) averaged over 5 replicates. Vertical and horizontal 
bars represent the standard deviations among replicates. Algorithm parameters used for APMC: a in {0.3,0.5,0.7} and Pacc„i„ in 
{0.01, 0.05,0.1,0.2}. Blue circles are used for Pacc,„„, = 0.01, orange triangles for Pacc„,i„ = 0.05, green squares for Pacc„j„ =0.1, and 
purple diamonds for Pace,,,;,, =0.2. PMC: red plain triangles for a sequence of tolerance levels from ej = 3 to 65 = 1 .4. SMC: grey plain 
square for (a,M) in {(0.9, 1),(0.99, 1)}, grey star for (a,M) = (0.9, 15) and a e target equal to 1.4. RSMC: brown plain diamond 
for a = 0.5 and a e target equal to 1.4. Results obtained with a standard rejection-based ABC algorithm are depicted with black plain 
circles, (b) Boxplot of the criterion "squared L2 distance times the number of simulations" for the different algorithms. APMC: for a 
in{0.3,0.5,0.7}andpacc,„„, = 0.01; SMC: for (a, M) in {(0.9, 1), (0.99, 1), (0.9, 15)} and a e target equal to 0.01; RSMC: for a = 0.5 
and a e target equal to 0.01; ABC: for a e target equal to 1.4: PMC: for a sequence of tolerance levels from ei = 3 to £5 = 1.4. (c) 
Criterion "squared L2 distance times the number of simulations" in the APMC algorithm for the different values of a and Pacc,„ii, ■ 
Each cell depicts the average of the criterion over the 5 performed replicates of the APMC. 
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5 Discussion 

The good performances of APMC should of course be confirmed on other examples. Nevertheless we argue 
that they are due to the main assets of our approach: 

- We choose an appropriate reweighting process instead of a MCMC kernel, which corrects the sampling 
bias without duplicating particles; 

- We define an easy to interpret stopping criterion that automatically defines the number of sequential 
steps. 

Therefore, we can have some confidence in the good performances of APMC on other examples. 

In the future, it would be interesting to evaluate this algorithm on models involving a larger number of 
parameters and/or multi-modal posterior distributions. Moreover, APMC could benefit from other improve- 
ments, in particular by performing a semi-automatic selection of informative summary statistics after the 
first ABC step (Joyce and Marjoram 2008; Fearnhead and Prangle 2011) and by using local regressions for 
post-processing the final posterior distribution (Beaumont et al. 2002; Blum and Frangois 2010). We did not 
perform such combinations in the present contribution, so that our algorithm is directly comparable with 
the three other sequential algorithms we looked at. However, they would be straightforward, because the 
different improvements concern different steps of the ABC procedure. 



Appendix A: Description of the algorithms 



Algorithm 1 Likelihood-free rejection sampler (ABC) 

Given N the number of particles 
for I = 1 to A' do 
repeat 

Generate 6* ~ ;r(6) 
Simulate x ~ f{x\0*) 
until p{S{x),S{y)) <e 

Set e, = e* 

end for 



Algorithm 2 Population Monte Carlo Approximate Bayesian Computation (PMC) 

Given W the number of particles and a decreasing sequence of tolerance level 6] > ... > 67-, 
For( = 1, 
for (■ = 1 N Ao 
repeat 

Simulate 6,*'' ~ ;r(e) and.r ~ /(x|e/'') 
unta p{S{x),S{y)) < El 

end for 

Take f7| as twice the weighted empirical variance of (6/'')i<,<a( 
for r = 2 to r do 
for ! = 1 to W do 
repeat 

Sample 9* from sj' with probabilities 

Generate s/'' 1 9* ~ .yK(e* , a^) and x - f{x\ s/'' ) 
until p(S(.v),5(v)) <e, 

Set »• ' ~ ^-r^ ' ^-7^ 

E5L,w<'-"arV(ar'(e,<"-(9j'-")) 



Take fT,^_^[ as twice the weighted empirical variance of (6/'')i<,<ai 



end for 

Take 
end for 

Where (p{x) = —=e 
•J2n 
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Algorithm 3 Sequential Monte Carlo Approximate Bayesian Computation Replenishment (RSMC) 

Given ei, ej, c,a€ [0, 1] andA'a = [aWj, 
for ! = 1 to do 
repeat 

Simulate 6, ~ 7t{9) and ;i: 9,) 
Pi = p{S(x),S(y)) 
until Pi < El 
end for 

Son {%, Pi) hy Pi 
Set Emax = Pn 
while Emax > Et do 

Remove the Na particles with largest p 

Set EneXT = PN-Na 

Set iacc — 

Compute the parameters of the proposal MCMC <?(■,■) with the N — Na particles, 
for = 1 to Na do 

Simulate dN-Na+j ~ {Si)l<i<N-Na 

for it = 1 i? do 

Generate 9* ^ q(9* ,6N-N^+j) lAx* 
Generate u < 'Wiff i 

Set dff-fla+j = 0* 
SetpAr-iv„+j=p(S(^*),S(3')) 

'acc ^ ^acc "t" 1 

end if 
end for 
end for 

I 

Set Pace — 



SetR= : 



RNa 

log(c) 



l0g(l -Pace) 

end while 
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Algorithm 4 Adaptive Sequential Monte Carlo Approximate Bayesian Computation (SMC) 

Given N,M,a e [0, 1], eo = °=, e and Nt, 

For? = 0, 

for (■ = 1 to A? do 

Simulate 0,'°' ~ 7t{e) 

for k=l M do 

Simulate z(;^j^/(.|e/°)) 

end for 

SetW'C"^- 
' N 

end for 

Wehave£5S((w/'''),eo) =iV where £'5S((W^.""),eo) = (i:£i(H'/''')2) ' 
Set ? = 1 

while e,_i > e do 

Determine £V resolving ESS{{wP),s,) = a£SS((H^.*'"''),e,_i) where « W;.''"'^ ''''''' '' et Ae,^, = 

{;(|p(5(x),5().))<£} 
if e, < e then 

£,j — £ 

end if 

if£S5((W^.'''),e,)<W7-then 
for ! = 1 to W do 

Simulate (e^;,"",^^'7,ll)) in (e*;-",^*;-;^,) with probabiUties , 1 < y<iV 

Setw(" = i 
end for 
end if 

for f = 1 to W do 
ifwj'' >Othen 

Generate 6* ~ ^:(e*le(|"'') 
for = 1 to M do 

Simulate JSf(.,j)-/(-|e*) 
end for 

Generate u < '^[o,i] 

i;f=iiA„,(X(,,));r(e*)*:,(e<'-"|e*) 

if M < 1 A ; — TT — -. — -. — 7— then 

Set (OflyXl^^i-M^) = 
Set = (^('^ 

end if 
end if 
end for 
end while 
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Algorithm 5 Adaptive Population Monte Carlo Approximate Bayesian Computation 

Given N, Na = [ciN] the number of particles to keep at each iteration among the N particles (a € [0, 1]) and Pacc„i„ the minimal 
acceptance rate, 
for r = 1 do 

for ! = 1 to W do 

Simulate e,*"' ~ K{e) and x - /(x|e/°') 

Setpf =p(5W,S(y)) 
SetH-,"« = l 
end for 

Let ei = Q (0) (a) the first a-quantile of p'"' where pC*' = (p/"' \ 

Let {(e/'\w«,pW)} = {{er,-r\pr)\pp < ^. ^ i < ' < 

Take crj as twice the weighted empirical variance of { (6/'', w-'')} !<;<?/„ 

Sit Pace = 1 

t<-t+l 
end for 

while Pace > Pacc„,i„ do 
for ; = A^a + 1 to W do 

,('-!) 

Pick from ef with probability p ,l<j<Na 
Generate e/'"'' | Of - ^(6? , af,_^s^ ) and ;c - /(;c| e/'"'') 



Setp<'-"=p(5W,S(y)) 



Set = ^(Q.''") 

Ej'^:(w<'-"/E£.wr'Vr-\9(cTr-Ue<'-"-er")) 

end for 

Set = iv=jv^ I*=JV„+l lp('-i)^£^_| 
Let e, = Gp(,-i) (a) where p('-') = {p*'"" 

Let {(e(",wW,p«)} = {{er\^-'\pj'-^\pr' <e.,l<i<N] 

Take of as twice the weighted empirical variance of {{Si'\wf^)}i<i<Na 
t-f-f+1 
end while 

Where Vh G [0, 1] and X = {;t:i, ...,;t„}, Gx(m) = inf{;c € Xlft W > u} and Fx(x) = ^ S=i l;ci<;c- 
Where (p{x) = t 
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Appendix B: Proof that the algorithm stops 

We know that there exists > such that £; Coo because, by constraction of the algorithm (e,) is a 
positive decreasing sequence and it is bounded by 0. 

For each e 0, we consider the distance {p{x,y)\d) as a random variable p(0). Let/p(e) be the probability 
density function of p(0). 

The probabiUty P[p(0) > £?] that the drawn distance associated to parameter 6 is higher than the current 
tolerance £( satisfies: 

¥\p{d)>e,] = l-V[{p(e)<et] 
= Itfp{e){x)dx 



We define: 



We have: 



^max = sup \ sup {/p(e) W} \ 

ee0 UeM+ J 



The N — Na particles are independent and identically distributed from 71;+ 1 the density defined by the 
algorithm, hence the probability P[pacc{t + 1) = 0] that no particle is accepted at step t + lis such that: 



N-Na 



V\pacc{t + 1) = 0] > {l-Vnuu{St - £00)) 
If 

^max < +°°> because £t — £00 — > 0, we have: 

nPacc{t + l)=0] 1 

We can conclude that Pacc{t) converges in probability towards if Fmax < +°°- This ensures that the algo- 
rithm stops, whatever the chosen value of Pacc„i„- 
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